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In order to use the Hess Structural Model to predict the need for certain cueing systems, George 
and Cardullo significantly expanded it by adding motion feedback to the model and incorporating 
models of the motion system dynamics, motion cueing algorithm and a vestibular system. This paper 
proposes a methodology to evaluate effectiveness of these innovations by performing a comparison 
analysis of the model performance with and without the expanded motion feedback. 

The proposed methodology is composed of two stages. The first stage involves fine-tuning 
parameters of the original Hess structural model in order to match the actual control behavior 
recorded during the experiments at NASA Visual Motion Simulator (VMS) facility. The parameter 
tuning procedure utilizes a new automated parameter identification technique, which was developed 
at the Man-Machine Systems Lab at SUNY Binghamton. In the second stage of the proposed 
methodology, an expanded motion feedback is added to the structural model. The resulting 
performance of the model is then compared to that of the original one. As proposed by Hess, metrics 
to evaluate the performance of the models include comparison against the crossover models 
standards imposed on the crossover frequency and phase margin of the overall man-machine system. 
Preliminary results indicate the advantage of having the model of the motion system and motion 
cueing incorporated into the model of the human operator. It is also demonstrated that the crossover 
frequency and the phase margin of the expanded model are well within the limits imposed by the 
crossover model. 


I. Introduction 

R ECENT research by Cardullo 1 et al (2006) significantly expanded the original Hess structural model 
by adding the vestibular paths as well as the models of the motion system and cueing algorithms. 
Figure 1 illustrates the expanded Hess model for both longitudinal and later channels. The highlighted 
boxes indicate the added motion feedback paths. Depending on which degree-of-freedom is being modeled, the 
expanded model utilizes models of otoliths or semicircular canals (SCC), which simulate motion perception in the 
translational and angular channels respectively. The motion feedback paths are then complemented by the motion 


1 Graduate student, Department of Mechanical Engineering, SUNY Binghamton, Binghamton NY, 13902, AIAA 
student member. 

2 Professor, Department of Mechanical Engineering, SUNY Binghamton, Binghamton NY, 13902, AIAA associate 
fellow. 

3 Research Associate, Man-Machine Systems Laboratory, SUNY Binghamton, Binghamton NY, 13902. 

4 Software Engineer, Unisys Corp, Hampton VA, 23666. 

1 

American Institute of Aeronautics and Astronautics 



platform dynamics and motion cueing algorithm. Platform dynamics is modeled by a second order transfer function 
(Eq.l): 








Figure 1. Expanded Hess model (Cardullo et al, 2006) with motion 
feedback for a) rotational channel, b) translational channel. 
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algorithm is represented by the classical high-pass filter (2): 
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With a natural frequency O) wo = 1 .02 ( ra ^/^ ) , and damping ratio ^ wo = 1.0 
The angular channel model utilizes the model of the SCC as proposed by Telban 2 et al (2005): 
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Where CO is the angular velocity and CO is the perceived angular velocity both in radians per second. The 
mathematical representation of the otoliths, used in the translational channel, was also developed by Telban et al 
(2005) and is given by: 


F p = (45 + .4) 

F ~ (.085 2 + 5.0165 + 1) ’ 


(4) 


Where F p is the perceived specific force and F is the actual specific force resulting from the translational 
motion. 

Parameters of all systems included in the enhanced Hess model are summarized in Table 1. Parameters 
inherited from the original version of the model were tuned manually so that the frequency response of the model 
complied with the requirements set by the crossover model in terms of the crossover frequency and phase margin. 


Table 1. Parameters of the expanded version of the Hess structural model 


Hess structural model parameters 

Control 
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K 

K e 

fa 
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K 2 



T o 
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i/ 
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1.0 

1.0 
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3.0 

2.0 

0.01 

- 

0.1 

0.707 

10.0 

mm 

1.0 

1.0 

20.0 

3.0 

2.0 

0.01 

0.1 

0.1 

0.707 

10.0 

Vestibular and Haptic rate feedback parameters 

Control Element 
Dynamics 

V 

see 

V 

^ Otolith 

K Haptic 

i/ 

As 

0.1 

0.1 

0.1 

1/ 

/V 

0.1 

0.1 

0.1 


The overall performance of the model with added sensory modalities was considered to be satisfactory: the 
phase margin with added feedback sensory data was improved, which was supported by the experimental data 
Cardullo et al (2006). 

One of the main goals of this investigation was to automate the process of tuning/identifying parameters of 
the Hess structural model with the emphasis on control behavior of a given individual. This paper addresses this 
issue by presenting the description of the Automated Parameter Identification algorithm (APID) as well as the 
description of the application of APID to aid evaluating the effectiveness of the motion/vestibular system feedback 
of the expanded Hess structural model. 
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II. Methodology 

The methodology used during this investigation to evaluate effectiveness of the motion/vestibular system 
feedback is based on evaluating the ability of the Hess structural model to simulate the behavior of an actual pilot 
with and without the motion/vestibular feedback. First, APID is applied to identify parameters of the original 
version of the Hess model, in order to get as close as possible to the actual control behavior of a given pilot. In the 
second step the motion system is added. The resulting change in performance is evaluated by comparing PSD data 
of the actual and simulated control signals. During the third step of the methodology, APID is applied to the 

expanded version of the Hess model, where gain K m is identified along with other parameters of the Hess structural 

model. The performance of the optimized expanded version of the model is then evaluated again by comparing the 
PSD data of the simulated control signal against that obtained in the previous step. 


III. Automated parameter identification algorithm 

The Automated Parameter Identification algorithm can be described as an automated optimization algorithm 
(Figure 2), which is searching for values of the identified parameters resulting in a maximum match between the 
actual and modeled operator control signals. The measure of “closeness” will be referred to as fitness of a particular 
combination of identified parameters of the model. The APID uses a genetic algorithm (GA) based optimization 
engine. This choice is dictated by several factors, which include: high rate of convergence, ability to deal with 
highly non-linear systems, and the ability to optimize a large number of parameters, which may or may not be 
related to each other. 



Figure 2. General architecture of the proposed APID 


A. Bit-Climbing Genetic Algorithm 

This section describes in detail one of the implementations of a classical genetic algorithm, namely the bit- 
climbing algorithm (BCGA), which was chosen to be used in APID. It was developed by Davis 3 (1991) and can be 
described as a modified hill-climbing algorithm with certain features adapted from the classical genetic algorithm. 
The BCGA uses binary representation of the string of estimated parameters and/or variables similar to the classical 
GA. However, unlike classical genetic algorithms, it requires only a single chromosome for operation. The name 
“bit-climber” is inspired by the fashion in which the algorithm manipulates individual bits of the chromosome. 
Davis discovered that BCGA converges to a solution from 3 to 23 times faster than a traditional GA, while 
maintaining an acceptable accuracy. Since computational efficiency is crucial for this research it was decided to use 
the BCGA instead of SGA or any other traditional GA. This was verified by implementing both BCGA and SGA in 
APID during preliminary research. 
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The bit-climbing algorithm begins with 
constructing the initial chromosome. Unlike 
the classical GA, the BCGA uses only a single 
chromosome, rather than a generation of 
chromosomes. On the other hand, the BCA 
chromosome is composed of identified 
parameters, put in a binary string, similarly to 
the classical GA. The initial chromosome is 
generated arbitrarily, i.e. there is no a-prior or 
“best guess” knowledge of the values of the 
identified parameters. Each identified 
parameters occupies 10 to 11 bits of the 
chromosome, depending on the desired 
accuracy. If, for example, the chromosome 
contains seven identified parameters and all 
parameters occupy 10 bits - it will consist of 
total of 70 bits. Once the initial chromosome is 
generated the iterative process begins by 
“flipping” (changing values from zero to one 
and vise-versa) bits of the chromosome in 
some arbitrary order. At each bit flip the 
fitness of the chromosome is calculated. Every 
time new chromosome produces a “better” 
fitness, the old chromosome is replaced with 
the new one and the bit flipping continues. 
When all bits are tested, proceed with the new 
iteration of flipping bits. The process continues 
a pre-allotted number of times or until no 
improvement in fitness is observed. Finally the 
“best” chromosome is saved. Figure 3 contains 
the flowchart of the bit-climbing algorithm. 

Every BCGA iteration of optimization is 
followed by the post-processing binary to 
decimal conversion, which may involve 
division by a factor of a 100. Hence, if a given 
parameter occupies 10 bits in the chromosome, 
it can take on values between 0 and 10.23 with 
the precision of up to a second digit after 
decimal point. In order to allow parameters to 
take on negative values, post-processing can 
include subtraction of a constant. Figure 4 
illustrates the post processing associated with 
identification of parameters of the Hess 
structural model. 


Figure 3. Bit-climbing genetic algorithm flowchart 


5 

American Institute of Aeronautics and Astronautics 












Figure 4. Pre- and post- processing associated with identification 
process within BCGA 

The BCGA optimization routine (the outer loop of the BCGA) is repeated for 10 iterations, thus ensuring the 
convergence of the algorithm to a solution. It was observed, however, that the BCGA algorithm usually converges 
within first three iterations. 

B. Metrics 

This section discusses the theory and application underlining the metrics used to evaluate fitness of the 
chromosome generated by the BCGA. In genetic programming, fitness score can be defined as an objective measure 
of how close the solution is to its maximum or minimum. The term solution used here means the optimum 
combination of identified parameters of a given model. Therefore, fitness of the chromosome is evaluated by 
comparing actual and simulated control signals in the power spectrum domain. The following paragraph discusses 
details of computing power spectral density of the available control signals. 

Power Spectral Density (PSD) of a control signal was chosen as an evaluation space for the control behavior 
since it allows quantifying the operator control activity in terms of the frequencies at which an individual operates as 
well as how much energy is spent at each frequency. 

The theory of signal processing, stipulates, that it is possible to estimate the PSD of a finite length signal data set 
with the help of the Fast Fourier Transform (FFT). The aim of spectral estimation is to express the distribution (over 

the frequency range) of the power in a signal. The spectral density of a fixed random process X n is related to the 
correlation sequence by the discrete Fourier transform and is given by: 


S.(®)=E (5) 

X=-oo 

Where R x is the autocorrelation sequence and CD is the specified frequency point. Eq. 5 can be rewritten as a 
function of the actual frequency f as follows: 


oo 

S,(J)= £ R,(x)e '■ 

*=-00 


( 6 ) 


Where f s is the sampling frequency. 
The PSD is then defined as: 
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(7) 


W) = 


s x (f) 


In order to reduce the error in spectral estimates, the method of averaging of the estimates over multiple 
segments of data is used. Such a method is also known as “windowing.” The original time history of the control 
signal is divided into a series of overlapping segments. Each segment is then weighted by the type of window used. 
In this analysis a Hamming window was used. This choice was dictated by the fact that it is superior to a strict 
rectangular window due to lesser side-lobe leakage. The length of each windowed segment was chosen to contain 
1024 data points. Overlapping regions contain only 32 data points. The Fourier transforms are then calculated for 
every weighted segment of data, which is further averaged over the resulting ensemble of PSD estimates. The 
resulting “smooth” spectrum estimate is then defined as: 


W)= 


f i A 


■ Un , i 


Zw). 


k = 1 


( 8 ) 


Where P x ( f ) is the power spectrum computed for each data segment and U is the correction factor associated 


with the energy loss and is equal to 0.612; n x is the number of overlapping segments and is given by: 


n 


X 



(9) 


Where T \ n - time duration of the spectral window; T' rec - length of the extended time-history that had been 

filled with the trim values; x frac - overlap fraction [ 0; l] (the lower value corresponds to “no overlap” and the 
upper value corresponding to 1 00% overlap) 

The fitness of a chromosome is computed by calculating the root mean square error (RMSE) between the two 
power spectrum estimates (Eq.10): 


E{f) = 


E(^ g (/)-w)) 2 


n 


( 10 ) 


In an attempt to quantify the fitness score space a benchmark score is introduced. The benchmark score is 
produced by executing the precision model, which is an extension of the McRuer crossover model. 

TV. Control task 

This section discusses an attempt to evaluate the effectiveness of the motion system feedback using actual data 

m T obtained at NASA Langley Research Center 

Table 2. Landing approach trim conditions. NASA LaRC VMS. (LaRC) yisua| Motion Simulator (V MS). 

Actual pilot control behavior was obtained 
during experiments which took place at 
NASA Langley VMS facility and simulated 
the landing approach of a Boeing 757 
aircraft under different control system delay 
conditions. One can find a detailed 
description of these experiments in the 
report by Guo 4 et al (2007). In this study 
only the basic landing approach is 


Altitude 

1300 ft BARO, 697 ft AGL 

Airspeed 

135 kts 

Heading angle 

180 deg 

Distance to runway 

2 nautical miles 

Flaps 

Full, Gear down 

EPR 

1.19 

Glide slope 

ON 

Localizer 

ON 
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considered, with no additional delay in the system. Table 2 contains trim conditions of the simulated aircraft. 

The control task that’s being investigated in this paper can be described as compensatory tracking of the roll 
rate of the aircraft in the presence of the lateral wind disturbance. The pilot’s task was to null the roll response of the 
aircraft due to the lateral wind. During experiments, the wind was simulated changing its direction clockwise as 
shown in figure 5.a. As the wind moves to the side of the aircraft it induces a sideslip angle of the aircraft. The latter 
causes the aircraft to roll due to dihedral effect. The side wind is also exerting some lateral force on the vertical tail 
of the aircraft thus contributing to the overall rolling motion of the aircraft (Figure 5.b). The lateral component of the 
wind can be recorded and used as a disturbance signal, which acts upon the aircraft. The human operator in this case 
perceives the rolling rate of the aircraft and attempts to null the effect of such disturbance by controlling the aircraft 
and counteracting the wind induced rolling motion. Figure 6 contains the block-diagram of the resulting man- 
machine system. 




Figure 5. 

a) Graphical representation of the wind change pattern during the 
simulated landing approach. 

b) Side wind acting on the vertical tail and causing the aircraft to roll 

8 

American Institute of Aeronautics and Astronautics 




Figure 6. Block diagram of a pilot-aircraft closed loop control system used to 
model the stabilization task in the roll channel 


In order to simulate the actual control behavior of a pilot during the lateral directional control of the aircraft it 
was decided to utilize a more sophisticated model of the aircraft, rather than a simple J/ or 2 dynamics. 

Roskam 5 (1979) offers a comprehensive description of the lateral directional dynamics of a wide fuselage jet 
aircraft. A complete transfer function relating aileron deflection 8 a to the roll angle (/> is given by Eq.l 1. 


(h A.s + B.s + C 

— = 1 t (11) 

8 a As 4 + Bs 3 + Cs 2 + Ds + E 

Where coefficients of the numerator and denominator are the functions of the lateral -directional dimensional 
stability derivatives of an aircraft and are given by the systems of equation (12) and (13) respectively: 

A = U(l-A l B 1 ) 

B = -Y f) (\-A,B i )-U(L p + N r + A,N p +B i L r ) 

C = U(L p N r - LN p ) + Yp (L +N r+ A x N p+ B X L ) - 

- Y p (Lp + NpA, + N Tf A,) + U(LpB l + Np + N Tfi )- 

’ - Y r (LpB l+ Np+N Tp ) ' (12) 

D = -Yp (. L p N r - L r N p )-Y p C L p N r -N p L r - N Tf L r ) - 

- g cos (0)(Lp + NpA x+ N Tp A) + (U — Y r )(LpN p -NpL p - N Tp L p ) 

[E = g cos (0)(LpN r - NpL r - N T0 L r ) 
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( 13 ) 


Aj=U(L Sa+ N s A,) 

Bt=U(N s L r -L s N r )-Y fi (L S ' +N S A 1 )+ 

{ +Y S '(L p + N p A + N T A x ) 

C,=-Y p {N s L r -L s N r )+\{L r N p +L r N Ti) -N r L p ) + 

{ +(U-Y r )(N p L Sa+ N Tfi L r -L p N Sa ) 

Where A 1 = and B 1 = . Table 3 contains expressions for the dimensional lateral-directional stability 

ZZ ZZ 

derivatives of an aircraft. 


Table 3. Lateral-directional dimensional stability derivatives (Roskam, 1979) 


Y - qSC » (ft/ ) 

fi ~ m ’ C Aec 2) 

qSC, 7 

L s = 4 , ( j/ 2 ) 

/ Vsec ' 

XX 

Y= qStC ” (//) 

p 2 mU ’ /sec 

N =- — ( V , ) 

^ / ’ Vsec 2 ' 1 

ZZ 

y = ?* C ,. (ft/ y 
r 2 mU ’ / sec 

qSb 2 C . 

N = P , (V ) 

p 21 U Vsec 7 

ZZ 

y - qSCys « (ft/ \ 

5a m ’ /sec 2 ' 

qSb 2 C . / 

2/ (7 Vsec 7 

ZZ 

qSbC, , / 

l b = ,p , (/ 2 ) 

p / v /sec 2 ^ 

XX 

qSbC 

N t = *>, (/ 2 ) 

v / /sec 

ZZ 

qSb 2 C, , / 

L - p ( V ) 

p 21 J7 ’ /sec 7 

qSbC , . 

/ V sec ' 

ZZ 

L _ <* Sb2 / ( 1/ , 

r 2I xx U ’ vsec ; 



1 T t2 

Where the parameter q is the dynamic pressure and is equal to q = — pv p 

Table 4 contains values of the dimensionless lateral-directional stability derivatives as well as geometrical 
parameters of an aircraft and flight conditions. 
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Table 4. Lateral-directional dimensionless stability derivatives (Roskam, 1979) and parameters of the aircraft 


Parameter 

Value 

Description (dimensions) 

s 

5500 

Wing area ( ft 1 ) 

b 

196 

Wing span {ft) 

U 

227 

Velocity (^/ c ) 

m 

564000 

Mass ( slug ) 

p 

0.002389 

Air density ( ) 

4 

13.7- 10 6 

Moments of inertia 

(slug- ft 2 ) 

4 

30.5- 10 6 

4 

43 . 1 - 1 0 6 

^ xz 

0.83- 10 6 


-1.08 

Lateral-directional 
dimensionless stability 
derivatives 

c „ 

0 

4, 

0 

C y Sa 

0 

c ’, 1 

-0.281 


-0.502 

4 

0.195 

c -. 

0.503 


0.184 

c.. 

-0.222 

C n 

n r 

-0.36 

c % 

0 

c -. 

0.0083 


For this research, however, the use of a reduced order aircraft transfer function (Roskam, 1979) relating roll rate 
(j) to aileron deflection 8 a was deemed to be sufficient (Eq.14). Such transfer function can be obtained as a result 

of “rolling” approximation of the lateral directional dynamics of the aircraft. The transfer function is derived under 
the assumption that there are no other lateral motion modes, such as dutch roll or spiral mode, being excited. 

i = — (14) 

s. S(s-L r ) 
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Substituting parameters from tables 3 and 4, equation 14 can be rewritten as: 


Aircraft.^ ^ 

a 5 a 4? + 0.0046) 


( 15 ) 


According to the block diagram shown on figure 6, knowledge of another model of the aircraft ( Aircraft w ), 

which would relate the roll rate of the aircraft to the side-wind, is required. Obtaining such model involved using the 
interactive Boeing 757 model provided by NASA Langley. Available aircraft (A/C) model allows simulating the 
behavior of the airplane in response to various pilot inputs and/or external disturbances, such as turbulence. The 
interface of the model allows choosing from a variety of predetermined control inputs, such as step, pulse doublet, 
sine wave etc. In addition to that, it is possible to feed any prerecorded, properly sampled external input as desired 
by the researcher. The modeled aircraft can be trimmed to virtually any flight conditions, including landing 
approach, cruise at low altitude, cruise at high altitude and so on. Recorded parameters match those recorded during 
actual simulation runs at the VMS research facility at NASA Langley research center. Availability of such 
interactive model of an aircraft allowed the derivation of the sought “side wind to roll rate” model by performing a 
standard parameter identification analysis based on the input-output data. The derivation process can be described as 
a four step process: 

Step 1. Trim aircraft to the initial conditions corresponding to the landing approach conditions used during 
experiments at NASA. 

Step 2. Feed the side wind external disturbance into the A/C model. Figure 7 graphically illustrates the side wind 
profile used during this derivation process. It is identical to that used during experiments at NASA. 



Figure 7. Lateral component of the wind acting upon the aircraft. 

Step 3. Record the roll rate response of the aircraft. 

Step 4. Using Auto Regressive with external (ARX) input model estimator, fit the available input/output 

data. 

With the procedure described above, it was possible to derive the following transfer function relating the 
lateral component of the wind acting upon the vertical tail of the aircraft to the roll rate of the aircraft (Eq. 16): 

</> 0.2554s 3 -\4.4s 1 + 212.6s - 25.72 

Aircraft w = — = — , (16) 

H,, s A + 84.375 s + 452.2 s 2 + 15 10s + 628.8 


12 

American Institute of Aeronautics and Astronautics 


V. Results 

Figure 8 contains PSD data of four signals: the blue solid line corresponds to the actual control activity of a pilot. 
The dotted black line represents the power spectral density of the pilot control behavior simulated by the original 
Hess structural model. Parameters of the model were optimized by the APID in an attempt to closely model the pilot 
control behavior. It can be seen, that the model performed quite poorly throughout the entire frequency range: the 
major power bin was underachieved substantially. 



Magenta colored dash-dotted 
line corresponds to performance of 
the Hess model with the motion 
feedback added into the system. 

Note, that the gain K m used in the 

model is equal to 150, which is 
much greater than the original 
value proposed by Cardullo et al 
(2006). If the original value of 0.1 
is used, then there is hardly any 
evidence of presence of the motion 
system feedback. At the initial 
stages of the investigation it has 

become clear that the gain K m 

must be increased in order to 
achieve some effect on the 
performance of the model. As one 
can see from the graph, the 


Figure 8. PSD of the actual control signal (blue solid) vs. simulated is u j • 

control signal by original Hess (black dotted), Expanded Hess m 

(magenta dot-dashed) and Expanded Hess tuned by APID (red noticeable change of the power 
dashed). spectral profile of the control 

signal: the power level around the 
main power bin of the actual 

control signal started to elevate. In the meantime, there has been no significant alterations of the power profile in the 


remainder of the spectrum. This observation suggests that inclusion of gain K m into a string of parameters 


identified by APID may improve the overall performance of the model. The fourth curve on the graph (red dashed 
curve), corresponds to the case when parameters of the Hess structural model were identified along with the gain of 
the motion system feedback. It is clear from the graph, that improvement was significant. The frequency (0.3 Hz) 
of the major bin of the actual pilot control signal was estimated very closely. Moreover, the power level was 
matched quite accurately as well. In the mean time, there has been no change in power spectral profile past 
0.8 Hz - it remained slightly elevated compared to the actual control signal. Table 5 contains numerical values of 
identified parameters of the Hess model for three cases considered. One can see only a slight improvement of the 
fitness score when the motion system feedback was added. The major improvement of the fitness score followed 

after APID had been applied to re-identify parameters of the Hess model including K m . 
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Table 5. Fitness scores produced by the original Hess structural model with and without the motion system 
feedback. 


Model 

configuration 

Identified parameters of the 
model 

Fitness 

K 


K 2 

K m 

Orig. 

16.2 

8 

2.08 

X 

0.0490 

Exp. 

16.2 

8 

2.08 

150 

0.0483 

Exp/ 

7 

3.35 

0.95 

976 

0.0263 


Application of the above described procedure to the control data from other pilots, who participated in the 
experiments demonstrate an average of 50% improvement in fitness score associated with introduction of the 
motion/vestibular system feedback along with re-tuning of the parameters with APID. 


VI. Conclusions 

This study was primarily aimed at evaluating the effectiveness of the motion/vestibular feedback developed 
by Cardullo et al (2006) for the Hess structural model of the human operator. The presented evaluation procedure 
utilizes a novel automated parameter identification technique (APID), which was designed specifically to aid in 
identification of parameters of a given model of the human operator. APID identifies parameters, which produce the 
closest approximation of actual operator control behavior. The uniqueness of the proposed APID is in using a 
genetic algorithm based optimization engine to achieve such a goal. The choice was dictated by the versatility and 
computational efficiency of genetic algorithms. There is also no obvious restriction on the number of identified 
parameters as well as there is no explicit requirement for them (parameters) to be related or influence each other. 

The proposed evaluation procedure was based on applying APID to tune parameters of the Hess model in 
order to model the control behavior of a certain individual. The actual data used during evaluation was obtained at 
NASA LaRC VMS research facility. Experiments were conducted with the active motion system. It was 
demonstrated that the Hess structural model is capable of close approximation of the actual control activity when it 
is enhanced by the motion/vestibular system feedback. Addition of the motion system feedback resulted on average 
a 50% increase in model performance. 

Further research in this direction may potentially contribute to answering the classical question whether one 
needs motion when simulating a certain control task. 
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